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ABSTRACT 

In this paper, we have numerically discussed the phenomenon of instabilities in a displacement process involving 
two immiscible liquids. The phenomenon is considered without involving the magnetic fluid. Numerical solution of 
governing non-linear partial differential equation for the phenomenon has been obtained by Finite element techniques. 
Finite element technique is a numerical method for finding an approximate solution of differential equation in finite region 
or domain. A Matlab code is developed to solve the problem and the numerical results are obtained at various time levels. 
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INTRODUCTION 

If a fluid contained in a porous medium is displaced by another fluid of lesser viscosity, then it is frequently 
observed that the displacing fluid has a strong tendency to protrude in form of fingers (instabilities) into more viscous 
fluid. This phenomenon is called fingering. 




OlSPLACEMtNl 



Figure 1: Fingering Process in a Porous Medium Oil Occupies 
Unshaded Area and Water Occupies Shades Area 

In Figure 1, fingering process has been shown between oil-water flows into a porous medium. In petroleum 
engineering, the fingering process is a well known phenomenon occurring in displacement of oil by water flooding that is a 
common oil recovery technique. 
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In the statistical treatment of fingering [1] only average cross-sectional area occupied by the fingers was observed 
while the size and shape of the individual fingers are neglected as in figure 2. Scheidegger and Johnson [21] discussed the 
statistical behavior in homogeneous porous media with capillary pressure. Verma [22] has examined the behavior of 
fingering in a displacement process through heterogeneous porous media. 



Displace meat 
Bjgctgd Liquid 



Native Liquid 



Figure 2: Schematic Representation of Fingers at Level V 

In the present paper, we have obtained a numerical solution of the problem by finite element techniques using 

Matlab. 

STATEMENT OF THE PROBLEM 

We consider that there is a uniform water injection into an oil saturated porous medium of homogeneous physical 
characteristics, such that the injecting water cuts through the oil formation and give rise to protuberance. This furnishes a 
well-developed fingers flow. Since the entire oil at the initial boundary (x=0) is displaced through a small distance due to 
the water injection. Therefore, we assume, further that complete water saturation exists at the initial boundary. 

MATHEMATICAL FORMULATION OF THE PROBLEM 

The seepage velocity of water (V w ) and oil (V 0 ) are given by Darcy's Law, 
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Where K is the permeability of the homogeneous medium, K w and K 0 are relative permeability of water and oil, 
which are functions of S w and S Q (S w and S Q are the saturation of water and oil) respectively, P w and P 0 are pressure of water 
and oil, 5 W and 5 D are constant kinematics viscosities, a is the inclination of the bed and g is acceleration due to gravity. 



Regarding the phase densities as constant, the equations of continuity of the two phases are: 
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Where, P is porosity of the medium. From the definition of phase saturation, it is evident that, S w + S 0 = 1 (5) 
The capillary pressure P c is defined as 

Pc = -PoS w (6) 

Pc=Pc-Pw (7) 

Where, (3 0 is a constant quantity. 

At this state, for definiteness of the mathematical analysis, we assume standard relationship due to Scheidegger 
and Johnson [21], between phase saturation and relative permeability as 



K w — S w 

K 0 = S 0 = 1 -S„ 



(8) 
(9) 



The equation of motion for saturation can be obtained by substituting the values of V w and V D from equations (1) 
and (2) into the equations (3) and (4) respectively, we get, 
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These are the general flow equation of the phase in homogeneous medium, when effects due to pressure 
discontinuity and gravity term in inclined porous medium are considered. 



Eliminating 
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from equations (10) and (7), we obtain 
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Combining equation (11) and (12) and using equation (5), we get 
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Integrating above eq. with respect to x, we have 
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Where, V is constant of integrating which can be evaluated later on. Simplification of (13) gives, 
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Using above equation in equation (12), we have 
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The value of pressure of oil (P 0 ) can be written as in [23] in the form 



p + p P -P - 1 
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(14) 
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Where, P is the mean pressure which is constant, therefore (15) implies 
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Using above equation in (13) we get, 
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Substituting the value of V from above equation in equation (14) we get, 
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Substituting the value of K w and P c from (6) and (8) we get, 



dS w P 0 K d 



dt 25,., dx 
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= 0 



(16) 



A set of suitable boundary conditions associated to problem (16) are 

S w (0,t) = 1 ; S w (x,0) = 0 ; S w (L,t) = 0 

Equation (16) is reduced to dimensionless form by setting 



X =x/L, T 



28, L P w 



, S (x,t)=S* (X,T) 



(17) 
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dS 



So that 
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With auxiliary conditions 

S w (0,T) = 1 ; S w (X,0) = 0 ; S W (1,T) = 0 

In equation (18) and (19) the asterisk are dropped for simplicity. 



(19) 



Equation (18) is desired nonlinear differential equation of motion for the flow of two immiscible liquids in 
homogeneous medium. 

A Matlab Code is prepared and executed with h =1/15 and k = 0.002223 for 225 time levels. The numerical value 
are shown by table. Curves indicating the behavior of saturation of water corresponding to various time periods. 

FINITE ELEMENT METHOD 

We attempt to solve the time dependent one -dimensional problem (18) by the application of finite element 
technique. We discuss the details of semi -discrete variational formulation of the problem. 

The domain of the problem, in present case, consists of all points between x = 0 and x = 1 Figure 3(a). 
This domain is divided into set of linear elements Figure 3(b). 
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Figure 3(a) 



h 

k ti-»j* I* H h h H 

■ m m w m m m 



m 7 [>U * DD H-1"fN~l 

Node Element 
Number Number 



Figure 3(b) 



Now, the variational form of given PDE (18) is 
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Choose an arbitrary linear element R (e) = [S/ e) , S 2 (e) ] & obtain interpolation function for R (e) using Lagrange 
interpolation Method such as 
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where N {e) =[N \ N 2 ] & 0 {e) =[S l S 2 f 

For linear element, Shape function is 
X - X l 
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Now, apply Variational Method to R (e) therefore equation (20) becomes 



R(e) 



+ 2S {e) — 



K dx , 



dT 



dX 



as, S {e \x) = N^ e) =0 ie)T N^ T 

■ _dN^ {e) _ {e) TdN^ T 

dX dX dX 

Therefore equation (23) becomes, 



and 



K dX y 



i{e) r dN ^ T avW ,( e ) 

dX dX 



(22) 



(23) 



/(s«)=± 



r dN [e)T 8N {e) ^ 
dX dX 



<f> (e) + 2\N {e >' N 



dT 



dX 



(24) 



For minimization, first differentiate equation (24) with respect to (j) 
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We use Gauss Legendre Quadrature Method to evaluate these integral. For this, first we transform co-ordinate X 
to a local coordinate z 3 when X = , z = -1; when X=S2,Z =1. 
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Therefore Interpolation function becomes iV x (z) = — (l— z) ', N 2 (z) = — (l + z) an d Jacobian matrix is 



J = — h Also element matrix transform to 
2 
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For A (e) , degree of polynomial p =2 then r = 2 and For B <e) , p = 1 then r =1. z Y and are corresponding gauss 
points and gauss weights with respect to 'r'. Then, the element matrix becomes, 
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ASSEMBLING OF ELEMENTS 

In deriving the element equations, we isolated a typical element (the eth element) from the mesh and formulated 
the variational problem and developed its finite element model. To obtain the finite element equations of the total problem, 
we must put the elements back into their original positions. The assembly of linear elements is carried out by imposing the 
following two conditions: 

• The continuity of primary variable requires 

s e =s e + l 
n 1 

• The balance of secondary variables at connecting nodes requires 
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S l if an external point surce of magnitudes \ is applied 



The inter-element continuity of primary variable can be imposed by simply renaming the variables of all elements 
connected to common node. For example, for a mesh of N linear finite element connected in series, we have 



For a uniform mesh of N elements, by equation (25) & equation (26), the assembled equation becomes, 

(27) 
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where, 
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Equation (27) represents the assembled equation. 

TIME APPROXIMATION 

We have obtained the finite element equation in the global form, which represent a system of simultaneous 
ordinary differential equation. We now introduced 8 family of approximations which approximates weighted average of a 
dependent variable of two consecutive time steps by linear interpolation of the values of the variable at the two time steps 



such as Sj =SS] +l +{\-S)S n j 



• S" +1 -S" 

The time derivatives 0 j are replaced by Forward finite difference formula such as, S j = 
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In view of (29) and (30), equation (27) written as, 
[A + 5 k(B((|) (n+1) )) ]()) <n+1) = [ A - (1-5) k(B(()) <n) )) ]§ (n) 
Where 5 =1/2 and n = 0,1,2, 

For a uniform mesh of N elements, by equation (28), the above equation takes the form, 



K{^ n + l) ) 



(n + 1) _ 



/») =F( /»)) 



(31) 



Equation (31) represents the global form of the problem. 

IMPOSING BOUNDARY CONDITIONS 

The boundary condition S w (0, T) = 1 states that Si (n+1) = 1 for all n > 0. Therefore we replace all the entries of 1 st 
row of matrix K by zero except the diagonal entry K u replaced by one and replace 1 st row of matrix F by 1 and boundary 
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condition S w (L, T) = 0 states that SN+i (n+1) = 0 for all n > 0. Therefore we replace all the entries of (N+l) th row of matrix K 
by zero except the diagonal entry K NN replaced by one and replace N+1)" 1 row of matrix F by 0. Thus, by applying 
boundary condition (19) to equation (31) and simplifying, we get 
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Thus, equation (32) is the resulting system of non linear algebraic equation. 

SOLUTION OF NON- ALGEBRAIC EQUATION 

In the previous section, we obtained the assembled equation which is nonlinear. The assembled nonlinear 
equations after imposing boundary conditions are given by equation (32). We seek an approximate solution by the 
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linearization which based on scheme 



Where, (j^' 1 ' denotes the solution of the n iteration. Thus, the coefficient Ky are evaluated using the solution 




from the previous iteration and the solution at the (n+1) iteration can be obtained by solving equation (33) using 



Gauss Elimination Method. At the beginning of the iteration (i.e. n=0), we assume the solution § from initial condition 
which requires to have Si (0) = S 2 (0, = = S N+ i (0) = 0. 

Graphical Representation 

A Matlab Code is prepared for 15 elements model and resulting equation (32) for N = 15 is solved by Gauss 
Elimination method. 

Saturation of injected liquid at time t = 0.1, t = 0.2, t = 0.3, t = 0.4 and t = 0.5 seconds are 
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Interpretation 

In all graphs, X-axis represents the saturation of injected liquid (S w ) of porous media of length one. Y-axis 
represents the time 't' in seconds. Solution obtained with h =1/15 and k = 0.002223 for 225 time levels. 

Initially, saturation of injected liquid is zero at each point on observed region. Also, there is full injected liquid 
saturation (i.e. S w = 1) at injected face x = 0 & there is no saturation of injected liquid at other end (x =1) irrespective of 
time. 

It is clear from graph that, for each value of T, Saturation Sw has a decreasing tendency along the space 
co-ordinate axis. Also, for each point of X, the Saturation increases as time increases but the rate at which it rises at each 
point in observed region slows down with increase in time. This shows that the stabilization of the fingers is truly possible 
with the assumptions made for capillary pressure and water saturation. 

CONCLUSIONS 

Engineers in several fields have to learn the mechanism of drainage and to apply to problems of water supply, 
land reclamation and stabilization of foundations and sub grade, and also to the fields of petroleum production and 
agriculture. 

Drainage in general is any provision for the removal of excess water. The common objective of land projects to 
prevent or eliminate either water logging or inundation or otherwise productive land. Drainage of projected land refers 
principally to the disposal of surplus natural water adversely affecting irrigation. Practically every area where irrigation has 
been carried on for time has been affected by high water table. Therefore provision for adequate drainage is an essential 
part of planning, construction and operation of an irrigation project. 

For agriculture purpose, the continued presence of water in excess of that needed for vegetation is harmful. 
Prolonged saturation of soil excludes air essentially for healthy plant growth and the soil becomes cold, sour and 
unproductive. Consequently unsaturated or irrigated soils is a necessary evil, so to this type of drainage where originally 
saturation conditions are existing up to the top. 
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